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ABSTRACT 


The aim of this thesis is to study the role of flood 
plains on flood peak subsidence and time of occurrence of peak .A tw< 
dimensional mathematical model of flood routing through a river with 
flood plains is formulated for this purpose. The model uses only 
continuity equation and the inertia terms are neglected . The 
two-dimensionality here does not refer to the unsteady flow equation! 
in two spatial dimensions , but rather to a physical situation in 
which river and the flood plains form a two-dimensional network in 
the horizontal plane . Two laws of exchange that govern the flow 
between adjacent cells are used in this model depending upon the 
natural topography . There are river type link and weir type link . j 
River type links are governed by the Manning's formula where as weir; 

""" I 

type links are governed by the weir discharge formulae . Bed slope j 

I 

is used in the Manning's equation for longitudinal direction where a4 

*> [ 

the water surface slope is used in place of bed slope for transverse | 

l 

direction. The results obtained from the present study are compared I 

i 

with the results obtained by Mahapatra (1990) using a one dimensional 
model . ; 

The various parameters that are considered in the present 
study are : 

j 

1. ratio of flood plain width to main channel width ( B p ) . | 

i 

2. ratio of peak flow to base flow ( Q r ) . ; 

3. ratio of flood plain roughness to main channel roughness ( N p )$ 

I 

(»i) : 



The rate of peak subsidense increases with an increase m 
value . This effect of Nr is more for higher values of . Rate 
peak subsidence increases with an increase in Br values . The rate 
flood peak subsidence decreases with an increase in B value 

r 

Difference in water levels in the transverse direction was found to 
be small when the water level was much above the flood plain level 
However significant difference in the water levels of flood 
plains and main river was found when the flow depth on the £1° 
plains is less than 10 


cm . 
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CHAPTER I 


INTRODUCTION 


L.l GENERAL 

Floods wreak havoc in many parts of the world and the 
situation is aggravated by urban and riverine developments .Much of 
:he trauma can be avoided by sound flood plain management , which 
rails for sophisticated planning tools . Among the most useful of 
these are mathematical models . The basic aim of such models is to 
simulate on a time scale the rise and fall of flood water 
Mathematical models can be used to 

(i) determine the extent of potential damage , 

(ii) design and operate hydraulic structures properly , 

(iii) predict the peak time , recession time and peak depths ,and, 

(iv) properly manage flood plains . 

Mathematical modelling of flood routing is basically the 
determination of hydrograph (stage or discharge) at any location in 
the channel if the hydrograph (stage or discharge ) is given at an 
upstream location . The two most important parameters in flood 
routing are (i) peak depth and (ii) time of occurrence of peak depth. 
Mathematical modelling of unsteady flow in rivets is much more 
complicated than in artificial channels because of possible 
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dimensional and 


inundation of flood plains . Flow may not be one - 
it is very difficult to quantify the topographical parameters of 
inundated areas . Many times it may not be possible to predict the 
extent of inundation itself , leave alone the other flood 
characteristics . Straight prismatic channels in which the flow may 
be considered to be strictly one - dimensional are seldom observed in 
nature .It is observed that over bank flows go their own way in many 
cases filling up the flood plains in a manner as dictated by the 
local topography . In some cases , over bank flow may never return to 
the main channel during the falling flood .Even if a channel with wel 
defined flood plains is considered , the flow pattern is not simple . 
It is seen that the interaction between the main channel and the 
flood plain flow affects the flood propagation significantly. 

Mathematical modelling of flood plain flow should consider 
the two - dimensional nature of the problem .Two - dimensionality 
here does not refer to the unsteady flow equations in two spatial 
dimensions , but rather to the physical situation in which channel and 
the storage cell form a two - dimensional network in the horizontal 
plane .As the flood wave propagates along the main river channel , 
the water from the main river cells will inundate the flood plain 
cells .Transfer of water can take place between adjacent flood plain 
cells also depending on the water levels .During the falling flood 
the reverse situation occurs i.e. the flow occurs from the flood 
plain cells towards the main channel, as the river stage lowers first. 
Most of the earlier models are one - dimensionally designed since 
two - dimensional models have historically been found to be either 
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too difficult or too expensive to construct and operate .However the 
lateral transfer of water due to transverse slopes influences the 
celerity and consequently the time of occurrence of peak depths 
(Cunge ,1975) .The flood is also attenuated because of this. 
Therefore , any accurate model should take into consideration these 
two - dimensional effects .In the present work ,a two - dimensional 
mathematical model is developed to study the role of flood plains on 
flood peak subsidence and time of occurrence of peak depth . 

1.2 REVIEW OF LITERATURE 

Mathematical modelling of unsteady flow using the 
complete Saint Venant equations is very common and is covered in 
most of the text books (Cunge et al 1980 ,Chaudhry 1993 ) .However, 
most of these methods have been applied to simple channels , and 
studies on flood wave movement through compound channels are only 
few . Verwey (1970) assumed that the overbank portion of flood plains 
serves only as a storage space and modified the usual one - 
dimensional Saint Venant equations accordingly . A distinction was 
made between the " storage width "used in the continuity equation and 
the " flow width " used in the dynamic equation . Yen (1978) also 
used a similar concept while studying the subsidence of peak flow in 
compound channels . He used a diffusion wave model for this purpose . 
Tingsanchali and Ackermann (1976) considered both the dynamic and 
storage effects of overbank flow in their one-dimensional model for 
flood routing in natural rivers .Effect of flood banks was considered 
through the momemtum correction factor , 0 and by using Lotter's 
method for evaluating compound channel resistance . It was observed 
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that higher peak stages and lower discharges are obtained when only 
storage effect of flood plains is considered . Later Tingsanchali 
and Lai (1988) using the same equations obtained semi -empirical 
exponential equations to describe the subsidence of flood peak 
discharge in channel during over bank flow periods . Tingsanchali and 
Ackermann (1976) and Tingsanchali and Lai (1988) used the governing 
equations in a non-conversation form and the Preissmann scheme for 
the numerical solution . Mahapatra (1990) and Rao et al . (1992) 

developed a similar one -dimensional model for flood routing through 
compound channels . However , they used the governing equations in 
conservation form which are simpler to integrate and a second - order 
accurate explicit finite - difference scheme for the numerical 
solution . Their numerical results indicated that the momemtum 
correction factor , (3 does not have significant effect on the flood 
peak subsidence and it can be taken to one . Mahapatra (1990) also 
concluded that the effect of inertial terms on flood peak subsidence 
is negligible and therefore , they can be neglected while routing the 
floods through compound channels . 

One - dimensional models for flood routing through natural 
channels with flood plains give erroneous results when the flood 
plains are vast . The interaction between the main channel flow and 
flood plain flow may result in lateral flow and transverse transfer 
of water . Mathematical models for flood plain flow should take into 
account these two - dimensional effects . Zanobetti et al . (1970) 
constructed and used a two - dimensional model for lower basin of the 
Mekong river . The model considered a system of ordinary first order 
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differential equations which represent water level as a function of 
time in some 3 00 cells in to which the whole area was divided . 
Exchange relations between cells were based either on simplified 
Saint Venant's dynamic equations or on the weir type exchange laws . 
Inertial terms were ignored in both of these exchange laws because of 
their insignificant effect . In two - dimensional models of flood 
plain flow the flood build up is relatively slow (except when the 
dykes break ) and the resistance terms (or friction slope ) prevail 
in the flow equations . It is then logical to drop the inertial terms 
from tidal equations because they are of small importance as compared 
to the water surface slope and friction slope . The system of 
ordinary first order differential equations was solved by an implicit 
finite - difference method . Later , Cunge et al . (1980) modified 
this model to take into account the rainfall and the evaporation and 
applied it to the river Senegal and to the Mono river basin . 
Hromadka et al . (1985) developed a simple two - dimensional dam-break 
model for flood plain study purposes . Approximate governing equation 
based on the diffusion wave approximation were solved using an 
irregular triangle element integrated finite -difference formulation. 
The model was applied to simulate a hypothetical dam-break over a 
flat plain . 


Akanbi and Katopodes (1987) developed a two - dimensional 
model for flood waves propagating on a dry bed . They solve the 
complete two - dimensional flow equation using a dissipative finite 
element technique . A deforming grid generation scheme was introduced 
to account for the effects of propagating or receding fronts . The 
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developed model was applied to a hypothetical problem involving a 
flood wave spreading radially on an impervious bed . However , no 
attempts were made to apply this model for simulating the spreading 
of flood waters on flood plains . Extension of their model for this 
purpose and to include the actual topographical information would 
involve some very elaborate procedures . Also , the computational 
requirements for the application of such a model are enormous. 
Leelerc et al. (1990) also presented a finite - element model for 
estuarian and river flows with moving boundaries . The model defines 
dry, partly dry and wet type of elements and adopts the governing 
equations for flow accordingly . The model was tested using a real 
life application to a tidal flow in the Maniconagan estuary in Canada 
Recently , Zhao et al . (1994) used a finite-volume method for 

developing a very versatile two - dimensional unsteady flow model for 
river basins . In their model , the river basin is discretized using 
a combination of unstructured triangular and quadrilateral grids . 
This model can deal with the wetting and drying processes for flood 
plain and wet land studies . This model was successfully applied to a 
portion of the Kissimmee river basin in Florida . However , this 
model requires small computational time steps , since it uses an 
explicit scheme for numerical integration . It should be noted here 
that the models developed by Akanbi and Katapodes (1982) , by 

Leelerc et al. (1990) and by Zhao et al . (1994) use the complete two 

- dimensional governing equations and require enormous computational 
power for their application . This level of sophistication may not be 
required for many engineering applications where the inertial terms 
in the governing equations do not have significant effect . In such 
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cases , the model described by Cunge et al . (1980) is very useful . 


1.3 SCOPE OF PRESENT STUDY 

Aim of the present study is 

(i) to present a two - dimensional mathematical model of a 
prismatic river with idealized flood plains , using an implicit finite 
difference method for the numerical solution . 

(ii) to study the role of various channel parameters on the 
flood wave propagation and the flood peak subsidence and. 

(iv) to compare the model with one - dimensional model using 

inertial effects developed by Mahapatra (1990) . 

1 . 4 CLOSURE 

The following chapters present the formulations and use of 
the two - dimensional mathematical model for flood routing through a 
river with flood plains . Chapter II presents the governing equations 
and the numerical methods to solve them . Chapter III describes the 
applications of the model to a prismatic compound channel and 
discusses the results . Conclusions and the recommendations for 
future work are presented in chapter IV . 



CHAPTER II 


GOVERNING EQUATIONS AND METHOD OF SOLUTION 

2.1 INTRODUCTION 

In this chapter , the governing equations for two 
dimensional storage routing of floods are presented . Numerical 
scheme for solving these equations is also described . 

Presentation in this chapter is organized as follows . In 
section (2.2) the various assumptions made during the derivation of 
the governing equations are stated . In section (2.3) -the governing 
equation is derived . In section (2.4) the laws of discharge between 
adjacent cells are discussed . Numerical scheme for the solution of 
the governing equations is presented in section (2.5) . 

2.2 ASSUMPTIONS 

The entire computational domain over which the flood 
routing computations are to be made is discretized, into several inter 
connected cells . The governing equation for the water level 
variation in each of the cells is derived based on the following 
assumptions . 

(i) Inertia terms are negligible due to gradual variations in 
the flow characteristics . 

(ii) The volume V of the water stored in cell is a function of 
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the height z s in the cell i.e. 

v . v (*,) 

(iii) The discharge Q" k between two adjacent cells i and k at 
the given time t^ - n.At is a function of the levels z" and z^ only 

i.e. Q n = Q fz n ,z n 

l,k * [ 1 ' kj 

(iv) The discharge between the two cells Q does not vary 

1 9 k * 

during the period of time At . 

(v) The slope , 8 of the channel bed is small so that 

sin 8 = tan 8 and cos 8 » 1 . 

(vi) The unsteady flow resistance co-efficient is assumed to be 
the same as for steady flow . 

(vii) The velocity distribution is uniform over the flow depth . 

(viii) The flow depth is measured at the center and it is uniforir 

throughout the cell . 


2.3 DERIVATION OF GOVERNING EQUATION 


In this section the continuity equation for a cell is 
derived . The derivation closely follows the derivation given by 
Cunge (1976) . 

A cell as shown in fig. 2.1 is considered . The water 
surface elevation in the cell above the datum is z t .A constant time 
step At is assumed . Then for any given time t R = n.At the water level 


in the cell i is z (t ) and the corresponding water surface ares 

I n 

(ABCD) is equal to A (t ) , the water surface area in the horizontal 

^ sin 


; [t i 
>i ”.*j 


at the new time 


plane of the cell . The water level is z 4 
level (n+1) .At and the water surface becomes A'B'C'D' . The surface 
area is equal to A (t ,) = A ft 1 + A A . The rainfall P (t) on 

n si I n+l J si ^ nj si 1 
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Fig 2.1 Continuity equation of a cell 
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the cell surface during the time At and the discharge Q z ,z 

i , k ^ 1 k j 

from the cells k adjacent to cell i being the reasons behind the 
increment of the water in the cell . The increase of water volume 
stored in the cell i during the time At may be considered in two ways 
from geometrical considerations 

2 (t + dt) 

AV i = J d2 t (2.1) 

Z i(t ) 


from discharge considerations 


AV 


I 


t + dt 
n 


P 1 (t) dt + 


i j 


t + dt 
n 


Q i.k ( V 2 k ) dt < 2 - 2 > 


where ( £ ) represents the summation of the volume coming in to cell 
i from all the cells adjacent to cell i in time At . Suppose that the 
surface A remains unchanged between the two levels and z «• Az s 


i .e . 


AA 


si 


A 


« 1 


8 } 


, then Eg .(2.1) reduces to 


AV i = A s i ( V Az i ------ (2.3) 

Comparing equation (2.1) and equation (2.2) one can write 

A gi (Zj) AZj > P s (T) . At + At ^ Q i tk ( z t ( T ) » z k ( * )) (2-4) 

k 

where , n.At s r * (n+1) At . 

If Az — > 0 , and At — > 0 , then equation (2.4) can be written 

in the differential form as 

d z ^ 

A .,-ar- - p , (t) + Y. Q ‘-“ ( V z *) <2 - 5> 

k 

This is the continuity equation for cell i .Similar equations can be 
written for all the other cells . The number of equations depends on 
the number of cells in the computational domain . For example , if 
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there are N cells then there will be N ordinary equations consisting 
of N , (z , z , - - - - - z ) unknowns . 

Id N 

2.4 LAWS OF DISCHARGE BETWEEN CELLS 

Two types of exchanges between cells are 
generally used .These are (1) River type link and (2) Weir type link. 
Mathematical representation of these links is described below . 


2.4.1 RIVER TYPE LINK 

These links represent the exchanges between two cells when 
there is no local obstacle to the flow and a mean resistance 
co-efficient for a given cross section of flow can be used . In the 
case of river type of links , Manning formula as given below is used. 

r.l/2 


l.k 


-1- . A . R 2 ' 3 
n 1 ,k 1 , k 


( 2 . 6 ) 


where 


l.k 


R. 


l.k 


n = Mannings roughness coefficient 

A. . = the area of flow cross section between cells 

i and k 

= the hydraulic radius of the flow cross 
section between cells i and k 
and S is the water surface slope 
In equation (2.6) , parameters A t k and R t k are functions 

of water level z in the flow section between cells i and k .This 

i , k 


enables one to write - 


4- v C - K ( Vk } 


(2.7) 


where , z = a z + (1 - a ) z and the function K - K ( z ) = 

l,k 1 k i.k 

K ( z , z ) is known as the conveyance factor of the flow section 

i k 
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between the cells i and k . a is known as the weighting co-efficient 
and is constant for a given pair of cells . It is generally 
calculated such that z^ is the result of a linear interpolation 
between levels z { and z^ . In the present study a is taken as 0.5 . 
The water surface slope in the equation (2.6) is approximated as 


S = 


z n+1 - z n+1 

k 1 

Ax 


( 2 . 8 ) 


where Ax is the fixed center to center distance between the cells i 
and k 

Now one can write a function 0 as 

, 2/3 


K (Z ) 

0(2 ) - 

l,k , — 

fix 


A R 


n V Ax 


(2.9) 


The discharge formula , whose sign will depend on the conventions 
adopted with regard to the direction of flow and may be written as 
follows - 

Qj k • sign (z fc - z { ) 0 (z ) fabs(z k - tT) (2.10) 

The function 0 ( z ) must be established a priori based on the 

* » K 

channel geometry . 


2.4.2 WEIR TYPE LINK 

In weir type links one uses the discharge formulae which are 
applicable to weirs . Two types of flow situations (Figs. 2.2 (i) and 

2.2 (ii) may occur for the cases considered in this study . The 
discharge equations for these are as follows . 

(i) if w <y (w+h) , the flow is considered to be free flow and 

the discharge equation is given as - 

Q cl . cdw . h 3/2 - - - - - (2.11) 

i,k 
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(i) During rising flood 



(ii) During falling flood 


Fig 2.2 Weir type link 
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where , cl= . 42 g . Ax 


( 2 . 12 ) 


(ii) w a — (w+h) , the flow is known as submerged flow an 


the discharge equation is given as - 


Q = c2 . cdo . w . h 1 

1 ,k 


(2.13) 


where c2 = 1.5 . cl .In the equations (2.11) and (2.13) cdw and cdo 
are coefficient of weir and coefficient of orifice and they are given 
by (Rouse 1936) 


cdo = 0.611 - 0.175 


f — 

{ h+w 


(2.14) 


cdw = 0.611 + 0.075 


(2.15) 


2.5 NUMERICAL SOLUTION OF GOVERNING EQUATIONS 

The governing equations are numerically solved to obtain 
the water level in any cell as a function of time .An implicit 
method is adopted in this study so that large values -of At could be 
used in the numerical solution . In this method , it is assumed that 
the discharge Q ( k | 2 ,, ( r ) > z k (t) j is an intermediate discharge between 

Q" and Q n * X . Now , the water level difference AZ in Eq.(2.4) 
cannot be solved explicitly because the right hand side terms contain 
the values of unknown functions at a time (n+1) At . Defining the 
intermediate value of discharge as 

e <c * (i - e > <2 - is > 

where , 0 s © £ 1 

Substitution of equation (2.16) in equation (2.4) leads to 

f > 

Az A * P At + At 6 V Q n *l + (1-0) Y Q?„ (2.17) 

i si i 1,1c i,k 

k k ' 

Discharge equations for the intercell exchange as given by equations 
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(2.6) , (2.11) and (2.13) are non-linear . Substitution of these 

equations in the equation (2.17) will result in a system of 
non-linear equations when written for the cells in the domain 
Taylor series expansion can be applied to linearize the equatio 
assuming that the variation in water level is small during the time 
interval At , making the method of solution easier . Taylor series 
approximation is applied as shown below 




- U - ©) Q" k « ° (2.18) 

k 

The unknown in Eq.(2.18) are Z n+1 and Z n+1 and Eq.(2.18) can be 

1 k ^ 

linearized as 



in which the subscript r indicates that the term in the brackets is 
evaluated at an iteration level r . AZ is equal to Z n+1 — Z n+1 . 

r+1 r 

which is the improvement in the assumed value of Z at iteration 
level r.In the equation (2.19) there are N unknown AZ values 
corresponding to N cells , thus making it impossible to solve 
explicitly . A system of linear algebraic equations for AZ t can be 
obtained by writing Eq.(2.19) for all cells in the domain . This 
system of linear algebraic equations is solved simultaneously for 
each time interval . A considerable amount of computer time would be 


16 




necessary if the system of linear equations (2.19) is solved 

simultaneously by any of the customary methods . However , the 

problem considered in this study has a special structure which can be 

exploited to reduce the computational time significantly . The 
variables Az t and Az^ actually involved in any single equation 

(2.19) concern only the cell itself and the four adjacent cells .This 
forms a five point molecule structure having many zero's in the 

matrix formed by the co - efficients of the variables Az 4 if a 
systematic numbering of cells is followed . such a system can be 
efficiently solved by the Strongly Implicit Procedure (Subroutine 
D03EBF ) available in the NAG Library . 

The above equations (2.19) are solved simultaneously for Az 
values . These values are added to the previous values of Z , such 
that the new values are Z = Z + Az . This procedure is repeated till 
the Az values are less than a tolerance specified by us . 
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CHAPTER 


TWO - DIMENSIONAL MATHEMATICAL MODEL FOR FLOOD 
ROUTING THROUGH A RIVER WITH FLOOD PLAINS 

3.1 GENERAL 

Unsteady flow in a river with flood plains is usually 
two-dimensional if the flow depth above the flood plains is shallow . 
In general these flows should be treated as two-dimensional flows for 
a precise analysis . By two dimensionality , we refer to the physical 
situation in which the cells form a two dimensional network . However 
the exchanges of water between the cells are purely one -dimensional . 
In this chapter , the two-dimensional implicit model is applied tc 
study the role of various parameters of flood plain characterizing 
the geometry and flow conditions during over bank flow periods . The 
results obtained from the present study are also compared with the 
results obtained from one - dimensional model developed by Mahapatrc 
(1990) . 

3.2 CHANNEL GEOMETRY 

A symmetric prismatic channel with flood plains as shown ir 
Fig. 3.1 is considered for study . This cross section is chooser 
primarily with a view to compare the results of the two - dimensional 
model with one - dimensional model studied by Mahapatra (1990) .The 
depth of the flood plains as measured from the bed of the river ,H | 
is taken as 5.0 m . There is no transverse bed slope in the mail 


18 




GO 

a 

*c5 

s 

"TU 

O 

O 

E' 

5 

• r- ( 

£ 


2 


<+H 

o 

a 

o 

•i—i 


> 

a) 


t/D 

00 

GO 

o 


a 


cn 

• ♦ 

bJD 

fL( 


19 


channel as well as in the flood plains . However , the longitudinal 
bed slope remains same for both the flood plains and the main channel 
The total length of reach is 90 km but results for only first 80 km 
have been studied in order to prevent the effect of downstream 
boundary c'oiuJit ion . Tho notations used in figure 3.1 represent 

B = width of the river 

mi 

B s = width of the river including flood plains 
= depth of flood plain bed from river bed 

3.2.1 IDENTIFICATION OF CELLS 

Although the results are represented only for the case 
of a symmetric prismatic channel , the model is actually developed 
for any size and shape of the flood plains and the river . The area 
vectors cross section of each cell and the center to center distance 
between adjacent cells need to be computed for the application of the 
two-dimensional model . For this purpose each cell is distinguished 
by a local node and the local node is in turn represented by its 
global nodes . In the present case there are 4 global nodes 
representing each cell placed at their four corners respectively . 
Once a particular cell i,j is considered then the local node of that 
cell can be easily prescribed . However , a problem arises with 
locating the global nodes of that particular cell . This problem is 
overcome by constructing a connectivity matrix , which connects the 
local nodes with the respective global nodes . For example , 
considering the figure (3.2) , there are 6 local nodes and 12 global 

nodes representing the 6 cells . The connectivity matrix of the 
figure (3.2) is given as 





1 2 6 5 ' 

2 3 7 6 

3 4 8 7 

5 6 10 9 

6 7 11 10 

7 8 12 11 

The connectivity matrix will contain as many rows as the 
number of cells present . The number of columns of the connectivity 
matrix is 4 . The coordinates of these global nodes are also 
specified as input to facilitate the computation of the area vectors 
of the intercell faces , plan area of each cell and center to center 
distance between adjacent cells . 

3.3 INITIAL CONDITIONS 

The depth of flow at the center of each cell is specified 
at t = 0 as the initial condition . A Steady uniform flow with water 
level slightly more than the level of the bottom of flood plains is 
assumed as the initial conditions . The exchanges of flow have been 
computed based on the river type link or the weir type link depending 
on the physical topography .River type of links are governed by 
Mannings equation where as the weir type of links are governed by the 
weir formula . The channel bottom slope is assumed to be constant 
along the computational domain and through out the duration of study. 

3.4 BOUNDARY CONDITIONS 
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3.4. 1 


UPSTREAM 


In order to compare the model under study with 
the one dimensional model developed by Mahapatra (1990) , a 
hypothetical input hydrograph having a log - pearson type III 
distribution with four parameters is selected for evaluating the 
effect of variations in the parameters . The inlow hydrograph is 
shown in figure (3.3) . The equation of the inflow hydrograph is - 

t 


Q(t) = Q + 
0 




- <t-t )/(t -t ) 
p g p 



p 


t - t 
<3 P 


where , 

Q b = base flow discharge (uniform discharge) 

Q peak = peak flow discharge 

t - time to reach the peak of inflow hydrograph 
t = time to center of gravity of inflow hydrograph 
t ~ time under consideration 


3.4.2 DOWNSTREAM 

It is assumed that the change in water level in downstream 
cells and the cells just preceding them is same . This is equivelant 
to prescribing a uniform flow condition at the downstream end and is 
only an approximation of what actually occurs in the field . This 
approximation procedure is many times resorted to because of 
non-availability of actual rating curves at the downstream end . The 
numerical results for the last 10 km of the computational domain are 
not considered in the analysis because of the above approximate 
procedure for the downstream boundary condition . 
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Discharge (m**3/s) 



Time (sec) 

Fig 3.3.1 Inflow hydrograph for Qr=4 



Fig 3.3.2 Inflow hydrograph for Qr = 2 
Fig 3.3 Inflow hydrographs 
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3.5 DEVELOPMENT OF THE MODEL 

The total area including the river and the flood plains is 
divided into number of cells . The task was made simpler by the 
assumption of well defined flood plains , otherwise the natura] 
topography oi. the area must be considered while breaking down the 
given area into cells . For the problem under consideration , there 
are four types of cells i.e. the river cells , next to river cells , 
intermediate cells and the side cells (next to side boundaries) 
Depending on the geographical topography each cell is linked with the 
adjacent cells either with a river type of link or with a weir type 
of link . There is an abrupt change in the bed levels between the 
river cells and the next to river cells .Therefore, they are connected 
with weir type of links and the remaining cells are all connected 
with river type of links . In the present case , the cell centers are 
the geometric centres , but depending on the cross section these will 
change . The direction of flow is allowed from centre to centre of 
adjacent cells and is determined with the help of area vectors . 

There are five flood plain cells on each side of the river 
in the transverse direction . The size of the river cells is 4500 m x 
50 m . The size of the flood plains in the transverse direction is 
dependent on the flood plain width considered in a particular run. The 
Manning's n and the water surface areas are not functions of water 
but remain constant .The side wall has the same Manning's n as that 
of flood plains . Bed slope is used along the longitudinal direction 
where as the water surface slope is used in the transverse direction 
while considering the river type link . 
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3.6 


COMPUTER PROGRAM 


Implicit method is used for solving the governing differential 
equation for water level variation in a cell . Newton- Raphson 
iteration method is used for solving the system of non-linear 
simultaneous algebraic equations . The flow chart is given below and 
the programme is enclosed in the appendix A . 

The main programme calculates the initial conditions , area 
vectors , sides of each cell , boundary conditions and solves the set 
of simultaneous equations by Newtons iterative method at each time 
step . As mentioned earlier , there are four types of cells depending 
on their position . They are river cells next to river cells , 
intermediate cells , and side cells . If one starts writing 
subroutines based on this classification , they will be lengthy . Any 
cell is connected to the adjacent cells through one of its four sides 
, hence it is simpler to write the subroutines based on this fact . 
For flow coming in or going out from top and bottom sides of a cell 
to its adjacent cells a seperate subroutine is written . Another 
subroutine is written for the flow coming in or going out from left 
and right sides through river link . Seperate subroutine is written 
for the weir type of links. Similar subroutines are written for 
calculating the elements of the jacobian . A seperate subroutine is 
written for converting these terms in to a matrix form . The elements 
of the jacobian matrix are later converted into the input form of the 
D03EBF routine in the main program . The computed stage hydrographs 
at 0 km ,40.5 km ,76.5 km measured center to center from the upstream 
cell are considered for the analysis . 
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FLOW CHART 


Read and write statements 
of channel geometry 


Calculate elements of 
connecivity matrix 


Subroutine 

Connect 


Calculate area vectors, 
sides , c/c distance to 
adjacent cells of each cell 


Subroutine 

Area 


Set initial conditions 


Set upstream boundary 
conditions 


set the values of water levels 
at next time interval 


Calculate the right hand side 
of the equation of the form 


[A] [X] 


[B ] 


Subroutines 

— Topbotm 
Ltrt 

- Weirlink 


Calculate the elements of the 
jacobian matrix 


Subroutines 

- Diff topbotm 

- Diff ltrt 
Dif fweirlink 

- Jacob 


Calculate dh i , dh - - 
using D03EBF 


- dh 















S' 
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3.7 


SELECTION OF INPUT DATA 


A rectangular symmetric , prismatic channel with well defined 
flood plains as shown in figure 3.1 is considered for the study .The 
length of the reach is taken as 90 km but the stage hydrographs are 
drawn up to a distance of 80 km in order to avoid the effect of 
downstream boundary condition . As shown in the figure 3.1 , the 

width of the main channel is taken as 50 m and the bank full depth as 
5.0 m . The Manning n for the main channel has been taken as 0.03 
where as the value of Mannings n for flood plains varied from 0.03 to 
0.06 . The longitudinal slope of the main channel and the flood plain 
are same and equal to 0.0005 . The initial uniform flow depth in the 
main channel is 5.5 and in the flood plains is taken as 0.5 m .The 
width of the flood plains is varied from 75 m to 375 m . The datum 
line is assumed to be horizontal and it is situated at 10 m below the 
bed of the upstream river cell . The total time of computation is 
taken as 1,50,000 sec i.e. 42 hrs and the time base of the inflow 
hydrograph is 70,Q00sec . The time to travel to peak depth at various 
downstream locations from the upstream end rate of subsidence depends 
on the channel parameters and the inflow hydrograph as 

T - f f x,h ,B , B ,rn,fn,g,Q ,Q ,t ,t ,T ] 
pl 1 ^ bins b p p g J 

i = f [ x, h , B , B , rn, fn, g, Q , Q ,t ,t ,T ] 

2 l bms bppgj 

where , i is the subsidence rate ,T is the time of travel of 

pi 

the flood peak from the initial station to the various downstream's 
stations , x is the distance along the channel , h is the initial 
uniform flow depth ,B s is the total width including both the flood 
plains and river ,B is the width of the main channel , rn and fn are 
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the Manning's roughness coefficients for river and flood plains 
respectively , g is the acceleration due to gravity , Q b is the base 
flow , Q is the peak flow , t and t are the time to reach peak 
and time to reach centre of gravity of the inflow hydrograph and T is 
the duration of the equivalent inflow hydrograph . In the above 
functions T is the time duration of an equivalent rectangle drawn 
having same volume as that of the inflow hydrograph as shown in 
Fig. 3. 3 . It can be written as 


T = 


pi 


r b 

I Q , (t > 


dt 


For comparison of the results with one dimensional model same 
formulae are used for normalising the peak depths and peak times . 
The parameter found for normalising the distance is k fe , 

A R 2/3 

” here - TT-TH- 

m 


Various non-dimensional quantities are defined below . 


The normalized distance 


x 


4 k T 

b 


The flood wave amplitude ratio = 


y - y 

J peak)x 2 in 1 tlal 

y - y 

J peak Inflow J Initial 


normalized time to peak depth is defined as T 


pi 


T , - 

pi 


time to peak. - time to peak. 

) x ^ ) Inflow station 


X 


f Q. 


+ 4 g hi 


where and A^ are calculated from the initial uniform depth 
in the river i.e. hir . 

The variation of parameters studied in this case are - 
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n = 

r 


fn 

rn 


Q - -J 
r Q. 


B = 

r 


B 

s 

TT 


These parameters are varied one by one and other parameters 


are kept constant . The results obtained in this study were compared 
with those obtained using one dimensional model . 


3.8 COMPUTATION 

The total channel reach is divided in to 2 0 cells in the 
longitudinal direction i.e. dx = 4 500 km and in to 11 cells in the 
transverse direct ion. The value of time step At was taken as 1000 sec. 
There were altogether eight sets of parameters studied in the 
present case , they are as follows - 



In all the above cases , the values of the following parameters 
remain unchanged unless and otherwise specified. 

B - 50 m g - 9.81 m/s 2 

m 

h = 5.0m hir = 5.5 m 

b 
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Ax = 4500 Im 
S = 0.0005 

O 

Q b = 565.5668 m 3 /s 
rn = 0.03 

Q = 565. 

i 


hif = 0.5 m 

t = 46,000 sec 
g 

t = 43,200 sec 

p 

tlast = 1,50,000 sec 
668 m 3 /s 


3.9 RESULTS AND DISCUSSION 

Static hydrogrnphs are plotted at 0 km , 40.5 km and 76.5 kt 

from upstream end in figures 3.4 - 3.11 for different paramete- 

combinations . Fig. 3.12 shows the water surface profiles along th< 

river at different times. These plots reveal that peak depth decrease; 
and time to reach peak depth increases along the river . Thi; 
qualitatively validates the two-dimensional model developed in th< 
present study . The results of the two - dimensional model were 
compared with that of one - dimensional model . 

The results of the numerical simulations to study th 

effect of Q , R and N on the peak depth are shown in the figure, 

r v r 

3.13 - 3.20 . The normalized distance is plotted on x axis and th 

flood wave amplitude ratio is plotted on y-axis in these figures . 
In figures 3.13 - 3.14 the effect of is seen on the flood pea! 

depth . The rate of flood peak subsidence for is observed to h 
more than that for B - 4 due to the availability of more storag< 

r 

in the flood plains . Similar trend is observed in the one 
dimensional model . In figures 3.15 - 3.17 the effect of N r i, 
observed . It in observed in all the three figures that higher th 
value of N higher is the rate of peak subsidence . Similar result 
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Fig 3.5 STAGE HYDROGRAPH AT DIFFERENT LOCATIONS ALONG 
THE RIVER 
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Fig 3.6 STAGE HYDROGRAPH AT DIFFERENT LOCATIONS ALONG 
THE REIVER 
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40.5 KM 



Fig 3.8 STAGE HYDROGRAPH AT DIFFERENT LOCATIONS ALONG 
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Fig 3.9 STAGE HY1 
THE RIVER 
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Fis 3.10 STAGE .HYDROGRAPH AT DIFFERENT LOCATIONS ALONG 
THE RIVER 
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Fig 3.12 WATER SURFACE PROFILES AT DIFFERENT TIMES(SECONDS) 
ALONG THE RIVER 




Fig 3.13 Effect of Br on Flood peak subsidence 




Qr = 2 



Fig 3.14 Effect of Br on Flood peak subsidence 
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NORMALIZED DISTANCE 

Fis 3.17 Effect of Nr on flood peak subsidence 




were obtained in one - dimensional model . Figure 3.18 shows the 
effect of different sets of parameters on flood peak subsidence .The 
maximum subsidence was observed for Q = 4 , B =16 and N = 2 and 

r r r 

the maximum rate of subsidence was observed for Q = 2 , B =16 and 

r r 

= 2 . The effect of can be observed from figure 3.18 . Higher 
values of result in higher flow depths above the flood plains . 

Higher flow depths result in more or less uniform flow over the 
entire cross section . Therefore , higher values of Q result in 

r 

lower rate of subsidence . Figures 3.19 - 3.20 show the variation of 
parameters for Q r = 4 and Q r = 2 respectively . In figure 3.19 , 
maximum rate of peak subsidence is observed for B = 16 and N = 2 

r r 

and in figure 3.20 , also the maximum rate of peak subsidence is 

observed for B = 16 and N = 2 . 

r r 

The effect of B^ on the time to reach peak depth at various 
locations along the river is shown in Figs. 3.21 and 3.22 . In these 
figures , the normalized distance is shown on the x-axis and the 
normalized time to reach peak depth is shown on the y-axis . It can 
be seen that higher the B^ value higher is the time to reach peak 
depth . Similar results are obtained from one -dimensional model. The 
storage effects of the flood plains are more for higher values of 
and therefore , it takes a long time for the peak to travel in the 
downstream direction . 


Stage hydrographs are plotted at 40.5 km from 
the upstream end in Figs. 3. 23 and 3.24 to see the effect of on 
peak depth subsidence .It is clearly seen that higher the value 
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Fig 3.19 Effect of different parameters for 
on Flood peak subsidence 
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Fig 3.20 Effect of various parameters on Flood pe 
subsidence for Qr=2 
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Fig 3.21 TIME TO REACH THE PEAK DEPTH ALONG THE RIVER 
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Fig 3.22 TIME TO REACH THE PEAK DEPTH ALONG THE RIVER 
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higher is the peak depth and also higher is the time to reach that 
peak depth . Figures 3.25 and 3.26 show the effect of B values on 

r 

the depth of flow . It is observed that with higher B values lower 
values of peak depths are obtained but the time to reach the peak 
depth is highei for higher values . Figures 3.27 and 3.28 show the 
effect of on the flow depth at 40.5 km station . It is natural to 
have higher values of peak depths with higher Q values and this is 

r 

observed in these figures . 


Water level was found to be almost equal in 
the flood plain cells and the river cells at any location along 
the river for any time in all the runs discussed earlier 
Computations were also made after decreasing the initial uniform flow 
depth . A maximum level difference of of 3 cm in the transvers 
direction was observed (Fig. 3.29 a ) when the initial flow depth 
over the flood plains was 1 cm . Water level difference in the 
transverse direction was 6 cm (Fig. 3.29 b) when the initial uniform 
flow depth over the flood plains was 1 mm . Further decrease of 
initial flow depth does not cause any more difference than 6 cm as 
shown in figures 3.30 (a) and 3.30 (b) . In figures 3.29 and 3.30 the 

transverse section is taken at 40.5 km from upstream end . It is also 
observed from the computations (Figs. 3.31 a and 3.31 b ) that as the 
:low depth in the flood plains become equal to 10 cm , the difference 
in the water levels along the transverse direction reduced to zero at 
my location along the river . Due to this reason , the transverse 
profiles in the figures 3.29 and 3.30 were drawn at 70,000 sec and 
lot at the peak stage . For the case of initial flow depth equal to 
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ig 3.26 Stage hydrograph at 40.5 km from upstream cell 

along the river showing the effect of Br on peak depth 
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Fig 3.27 stage hydrograph at 40.5 km from upstream cell 
showing the effect of Qr on peak depth 




Depth above flood plain bed (metres) Depth above fi 00 d plain bed (metres) 



(a) with initial unifom flow depth of 1cm above flood plain 



(b)wilh initial uniform flow depth of 1mm above flood plain 

Fig 3.29 TRANSVERSE PROFILE AT 40,5 KM 

FROM UPSTREAM CELL Qr=4 Br=16 Nr=2 











Distance (metres) from left to right 

Fig 3.31 (b) Transverse profile at 40.5 km from the upstream cells along the river 
at various times for initial flow depth of 5.0001 metres . 




be 5.5 m , the water level difference in the transverse direction was 
observed to be almost equal to zero . The effect of initial flow 
depth on the flood peak subsidence is shown in Fig. 3.32 . Lower 
initial flow depths introduce significant two-dimensional effects and 
consequently increase the rate of flood peak subsidence. The two - 
dimensional effect on the peak depth translation in the downstream 
direction is shown in Fig. 3.33 . 
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Fig 3.3 '2 Effect of initial flow depth on 
Flood Peak Subsidence 





CHAPTER IV 


CONCLUSIONS 


4.1 CONCLUSIONS 

Various conclusions made after the study of the results are 
as follows - 

(i) Results obtained from two - dimensional implicit model for 

flood routing without considering the inertia terms slightly differ 
from the one - dimensional model results developed by Mahapatra .K. . 

(ii) Almost similar results were obtained with one 

dimensional model for peak subsidence with Q r = 4 and = 2 and 

B = 4,16 . 

r 

(iii) Flood peak subsidence was observed for all sets of 
parameters .Time to reach peak depth increases with an increase in B r 
values . 

(iv) The observed difference between the flood plains and the 
river cells was not much . Runs were made with initial uniform flow 
depth up to 5.00001 m and it was observed that as soon as the depth 
in the flood plain cells becomes 10 cm high , then the water levels 
are almost found to be equal in flood plain cells and river cells . 

(v) Seeing the variation of different parameters on the flood 
peak subsidence , it is concluded that the maximum rate of flood peak 
subsidence is observed with B = 16 and N = 2 . 

r r 
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4.2 RECOMMENDATIONS 

(i) The same case may be studied by taking the initial flow depth 
below the flood bank .i.e. the flow is only in the river . 

(ii) In the development of the model the physical features of the 
flood plains should be considered such as roads , dikes etc. 
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APPENDIX - A 


TUO DIMENSIONAL MATHEMAT I CL A MODEL FOR FLOOD ROUTING 

L-NUMBER OF CELLS IN LONGITUDINAL DIRECTION 

N-NUMBFR OF CEILS IN THF TRANSVERSE DIRECTION 

JN , NP *• NUMBK R OF CELLS IN THE COHI’Ul A I I OMAL DOMAIN 

T-TIME UNlirR CONSIDERATION 

DT" INCREMENT IN TIME INTERVAL 

HD-DEPTH OF BED OF UPSTREAM RIVER CELL FROM DATUM \ 

HM»HE I GUT OF THE VIETTED SURFACE IN RIVER 

HB"BEPTH OF FLOOD BANK FROM RIVER BED 

DZCHK«M IN IMUM ALLOWABLE CHANGE IN WATER LEVEL 

DZMAX=MAX I MUM CHANGE IN WATER LEVEL 

TLAST=END OF THE COMPUTAT I ONAL TIME 

TPEAK=TIME TO REACH PEAK DEPTH OF INFLOW HYDROGRAPH 

TGRA*TIME OF CENTRE OF GRAVITY OF INFLOW HYDROGRAPH 

RN=MANMING'5 ROUGHNESS COEFFICIENT OF RIVER 

FN“MANN U4G ' S ROUGHNESS COEFFICIENT OF FLOOD PLAIN 

HIR-INITIAL UNIFORM FLOW DEPTH OF RIVER 

HIF»INIT I AL 'UNIFORM FLOW DEPTH OF FLOOD PLAIN 

QD®R AT 1 0 OF PEAK DISCHARGE TO BASE FLOW DISCHARGE 

XF 1 U , XFEU*X- COORD I NATES OF FLOOD PLAIN OT4 UPSTREAM SIDE 

XF1D,XFBD=X-COORDI NATES OF FLOOD PLAIN ON DOWNSTREAM SIDE 

YF 1 U . YF2U = Y-COORD I NA T ES OF FLOOD PLAIN ON UPSTREAM SIDE 

YF l 0 i YP St'* Y-C OORD I NATES OF FLOOD PLAIN ON DOWNSTREAM SIDE 

XRTU.xnaU’X-COOPOtNATF.S OF RIVFP ON UFSTREAM SIDE 

XRID, Xfian-X-COOROINArr:S or-- rtvi-r on DOWNSTREAM SIDE 

YR 1 U , YR Y - COORD I NA T ES OF RIVER ON UFSTREAM SIDE 

YR1D, YE2D“Y-C0DR0IMATES OF RIVER ON DOWNSTREAM SIDE 

HI “DEPTH OF WATER LEVEL FROM BED AT KNOWN TIME 

HEADER T H OF WATER LEVEL FROM EED AT NEW TIME 

ZT=DEPTH OF WAFER LEVEL FROM DATUM AT KNOWN TIME 

ZE=DEPTH OF WATER LEVEL FROM DATUM AT NEW TIME 

aa*souRCE terth net discharge) 

A 1 “ELEMENT C OF THE JACOBIAN MATRIX 

BT.8L.HP ER *LENG THE OF THE TOP , I.FF T . b'O 1 TOM AND RIGHT SIDES OF EACH CELL 
DXT , DXt. . UXR , DXP»C/C DISTANCE OF CEILS ON TOP , LEF T, RIGHT AND BOTTOM OF 
ANY PARTICULAR CELL 

CB, CL . OR . CT“ARE THE CORRESPONDING AREA VECTORS TO ANT CELL FROM BOTTOM, 
l EFT, RIGHT AND TOP CELLS 

URKOP 1 , URKSPE , URKSP3“UQKK SPACE PROVIDED FOR NAG ROUTINE 

RESI05, CHNGSi APAPAH, CONRES, CONCHN* GOVEPM 1 HE COVERGENCE OF NAG ROUTINE 

AS GIVEN IN NAG MANUAL 
A, B , C , D . E“t NPU T S TO ! HE NAG ROUTINE 

X = CHANSE IN THE U A ’’ E R LEVELS AFTER EVERY ITERATIONS 
IMPLICIT DOUBLE PRC- C I S I ON ( A -H , 0- Z ) 

EXTERNAL D03EBF 
P AR AME T ER t MF 0 c 09 ) 

D I MENS I ON H 1 (1 0 0 , 1 1 ). He (TOO, IT ),Z1(TOO,11),ZEtlOO,1T),A1tNP,NP), 
j QQ< t 00 , I 1 ) , BH 1500 ) , BL M 500 ) , BB( I 50 0 > , BR U 500 ) , CR < 1 500 ) , 

1 CB( TSOO),CL( 1500) ,DXI_( 1500) . DXR ( 1500) ,DXT(1500) , DXB ( 1 5 0 0 ) .. 

1 , URKSPT ( TOO, T 1 ) , WRK5P2 ( TOO, 1,1 ) , URKSP3( 100, J I ) , RES IDS ( 250 ) 

* , CHNGSt 250 ) ,CT « 1 500 ), AC( 1500 ) , 

* A{ioo,in,Bnoo,ti),ctioo,ii),o(ioo,ii),E<ioo,ii),xnoo,ii) 

OPEN ( UNI T“ 1 , F IL.E“ ‘ IND ' ) 

OPEN ( UNI T=2 , F ILE= ' OUT DT ’ ) 

READ (1 , *■ ) L , N , JN . DT . HD . HB , DZCHK 

READ! T , ♦ ) TLAST, Tpr AIT, TGRA.HUI ,HUP,HUB 

READM.** RN.FN.SO H IR , HIF . OD . AFARAM , CONRES , CONCHN 

RE AD ( 1 i * 1 XF 1 U • XR 1 U , XR HU , XPe'U , XT I D , XR 1 D , XR 2D , XF 2D 

READ! 1 , * ) YFTU.YPTU, YR2U, YF2U, T F 1 D . Y R 1 D, YR£D, YF2D 

CALL AREA (N , XF T U ■ XP 1 U , XR2U, XF£U , XF 1 D , XR 1 D , XR2D , XF2D , 

1 YF1U.YP1U . Y R 2 U • Y F ? U . (F 1 D . YR T D , TRPD . YFED , 

I BT , BL , BB, BR , DX T , D>’l , DXB , OXR , CT , CL , CB , C R , AC ) 

DX»XF1 O/M 

BRA“( YFSU-YFTU)/t YR2U-YRTU) 


NRA*rN/PN 

WRITE(2, ' ) ' HR = ' .»RA. ' 


QU,' BR “ * , BRA , DX * *,DX 



T*> 0 . 0 

DOIU'1,1 

DO 1 1 1*1 ,N 

IF< J . EG . 6 J THEN 

HI ( I , J)»H1R 

Z 1 ( I , J ) =H 1 ( I ,J)+HD-( 1-1 ) * DX « SO 

ELSE 1 

HH J, Jl'HIK 

ZUI,JI«H1U.J)+HOUIB-.l 1-1 ) *DXt SO 
END 1 F 

1 1 CONTINUE 

IF1HIR LT .MB) THEN 

HM-HI R 

ELSE 

HM=HB 

ENDIF 

OS* t BEK U *HIF ) + * t 5 . /3 . J+SORM SO )/( FN+ <BB( U+HIF )+* 1 2. /3. 1) 
QI=BB<N+I )+HIF + ! M5./3. ) •» SOR 7 1 SO ) /F N 
QN=BB M . 0 + N + 1 )»HIF*+(5,/3. ) * SOFT T 1 50 1 /FH 
QF 1 * ( BB ( 5 . 0 + N + t ) +HI R ) *+( 5 . /3 ) *SORT ( SO) 

QF£*= ( RN * ( BB ( 5 . O + N+l 1 < 3 . 0 ♦HR ) ( 2 /3 . ) J 

QF = QF I /'OF 2 

5 T=T+DT 

DO 1 2 J = 1 , L 
D0 1 £1=1 ,N 
H2 U , J ) =M 1 (I , J 1 

1 2 CONTINUE 
TOIF**TGRA-TPEAK 

QM*QF+t OD-1 ) TQF + EXP ( ( TPEAK-T )/T DIF )* { T/TPEAK J + * ( TPEAK/TD I F ) 

T DX=XF1D/» 

D06 J - 1 , L 

DOG I = 1 , N 

I F ( J . EO . 6 ) THEN 

HB=0 . 0 

ELSE 

MB "5 . 0 

ENDIF 

6 Z2( I , J ) = H£( I . J M-HD + HB- ( I -1 H DX+SO 
D013J- I ,t 

D01 31=1 . N- 1 
IE=< J-1 J + N + I 
IB»1E+1 -J 

0 C=AC ( IE (HSU , J)-H1 ( I , J) )/DT 

CALL TOPBOTMU, J. IE,FN,RN.Z1 ,Z2,H2,HB,BT 
1 , BB,CT,CB,DXT,DXB,QS,QI , ON , QM , QT , QB ) 

IF( J.EQ . 1 ) THEN 

CALL LTRTU, J, J + 1, IE.FN.Z1 , Z2 , H2, BR, CR, DXR , Q ) 

F=QC-QI-QB-Q 

ENDIF 

IF(J.GT 1 . AND . J . L T 5 OR.J.GT 7 AND . J . LT . 1 11 THEN 
CALL LTRT (I , J, J + 1 , IE ,FN,Z1 , Z2 , H2, BR , CR , DXR, OR) 

CALL LTRT ( I , J , J - 1 , IE . FN , Z 1 , ZP , H? , BL , CL , DXL.OL ) 

F = QC.-Q 1 -QB-QR - QL ’ 

ENDIF 

I F < J EO . 5 ) T HEN 

CALL LTRT ( I , J, J-1 , I E » FN , 7 1 , Z2 , H2 , BL , CL , DXL , CJL ) 

CALL UEIRLINFU, 1+1 , I , IE , BR , CR . H 1 , H2 , HB , QR ) 

F=QC-QT-OB+OR-QL 

ENDIF 

IF( J.EQ. G) THEM 

CALL UEIRLIHKt I , J , J - 1 , I E , BL , CL , HI , H2 , HB , OL 1 
CALL UEIRUN!!( I . J , J+ 1 , IE, BR , CR , HI , H£ , HB , QR ) 

F=QC~QT-QB-OL -OP 
ENDIF 

1 F ( J . E 0 7 1 I HEN 

CALL LTRTU, J, J+1 , IE ,RJ,Z 1 . Z2, M2, BR, CR,DXR ,0R) 


72 



CALL UEIRLINKt I . J-1 , J , IE,BL , CL, HI , HE, HB, OL) 

F«QC“Or-aB-QR+OL 

ENDIF 

IFCJ.EQ.il ) THEN 

CALL LTRT( I , J ,J-1 , 1 E , FN , Z l , Z2 , H2 , BL , CL , DXL ( Q ) 

FaQC-QT-OB-Q 
ENDIF 
QQCI, J)*F 
CONTINUE 

CALL JACOBIN, JN, FN ■ R N , Z I » Z 2 , H 1 , H 2 , DT , AC , DX T , DXL , DXB , DXR , 
BT,8L,BB, BR, Cl , CL, CB . CR.HB, A! ) 

DZMAX*0 0 

DO J-1 , L 

DO I ■= 1 , N- 1 

J1 = C J-1 )+ CN-1 ) +•] 

CU,J) ! 'A1(J1,J1) 

IFCJ.GT 1)1 HEN 

A( I , J ) = A 1 ( J 1 ,J1-N+1 ) 

ENDIF 

IFd.CT 1 ) THEN 

B ( I , J 1 R A 1 t J 1 , J t - 1 J 

ENDIF 

IF (I .LT.N-1 ) THEN 
DC I, J ) * A 1 ( J1 . J1 + 1 ) 

ENDIF 

I F C J LI 1 I )T HEN 
E 1 I , J )=-Al ( J 1 . J 1 +N- 1 ) 

ENDIF 

ENDDO 

ENUDO 

DO 

DO J" 1 , L 
X ( I , J ) = 0 0 
ENDDO 
ENDDO * 

CALL DO JEFF I 17, 11, 1 0 0 , A , B , C , D , E , 00, X, AFAR AM ,250,0, ITU5ED ,1,1,1 
, C0NRE5 , CONCHN. RES IDS . CHNG5, IJRKSP1 , WRKSPE , WRKSP3 , 0 ) 

DO 1 5 1 = 1 , N - 1 
DO 1 5J = 1 , L 

IF(APS(X(I,J1) I. T . DZ11AX ) T HEN 

DZHAX=DZMAX 

ELSE 

DZMAX=ABS( XC I , 3 ) ) 

ENDIF 
CONTINUE 
DO 16 J = 1,L 
DO 161 = 1 ,N-1 

HSU , J ) «H2( I , J > -XU , J ) 

D01 7J“1 , L 

H2 ( N , J ) =H2 < N- 1 , J) 

UR I T E (*■,+■ ) ' DZMAX* ’ , DZMAX, * T»',T, ' ITUSED «',ITUSED 

IFCDZMAX. GT . PZCHK ) GO 107 
UR I T E t ? , * ) " 7 ~ , T 

WRITE CB, 1 9 ) 1 (H£f I . J ) , J = 1 , L) , 1=1 ,li) 

FORMA K5F7 5,2X,6F7.5) 

D0 1 8 J = 1 ,1. 

D01 81=1 ,N 

HI U , J ) =H2tI , J ) 

DO M J = 1 , L 
D0 1 A I -1 , N 
IFCJ.EO C' 1 THEN 

Z1 U , J) --Ml U , J 1 +HD * t I -1 JUUtJO 
ELSE 

MUX **50 

ENDIF 

CONTINUE 



IF< T. LT . TLAS7 ) G0T05 

STOP 

END 

C ELEMENTS OF CONNECTIVITY MATRIX CONNECTING 

C GLOBAL MOOES TO LOCAL NODES, 

SUBROUTINE CONNE CT < N , C ) 

IMPLICIT DOUBLE PRECISION !A-H,0-Z) ' 

DIMENSION C( 1 BOO . 1 BOO) 

D03J=1 , 4 

IF! J. FS. 1 ) T HEN 

D04 1= 1 , N 

4 C t I , J 1 = I 

DOS I *• N M ,2>N 

5 C ( I , J ) = I ♦ I 
DOGI=2+N+1 ,3«N 

6 C(I,JJ-I+2 
DOT I - S'* N + 1 , 4 *N 

T CC I , J)=I+3 

DOB I = 4+ N+ I , 5'f N 
B C ( I , J ) ~ I +4 

D09 1=5TN+1 , G+N 
9 C( I,J > = I -* 5 

DO 181=6+ NH , 7 * N 

18 Ct I , J )=I»6 
D022 I =7+N+ 1 , •?. 4 N 

SS CCI,J)*I+T 

D0£3I=S+N+1 , 5 *N 

53 CII,J)'l 4 8 
DO£4I=?+N4 1 , I 0 * N 

54 Ctl,J)=l+9 

DOS51 -1 0 »NM , 11 IN 

55 C ( I , J ) - I M 0 
END IF 

IF! J .FO 1 SIT HEN 
DO 191 = 1 . J 1 < N 

19 C ( I , J ) = Ct I , J-1 ' + 1 
END IF 

I F ! J EO . 3 ) THEN 
D 0 2 0 I = 1 , 1 1» N 

SO C( I, J) = C< I, J-S1+N+S 

ENDIF 

I F ( J EQ . 4 ) THEN 
D021I=1,U+N 

21 C(I,J) = C<I,J--1)-1 

ENDIF 

3 CONTINUE 

RETURN 
END 

C SUBROUTINE TO CALCULATE THE AREA OF 

C EACH CELL AND GIVE THE CORRESPONDING AREA VECTORS. 

SUBROUT I ME ARE A ( N , XF 1 U , XR 1 U , XR2U , XFSU , XF 1 D , XR 1 D , XRSD , 

1 XFSD . Y F 1 U , YR 1U . YR2U , YF2U, YF1 I) . YR1 D , YRSD , YF2D , BT , BL , BB , 

1 BR , DXT , DXt. , DXB , DXR , Cl , CL , CB , CR , AC ) 

IMPLICIT DOUFLE F R EC I SI ON ( A-H , O-Z ) 

DIMENSION XU 50 0) , Y( 1500) ,EXf 50 ) ,EY( SO) , CM 500, 1500) ,AC< 1500) 

1 , BTU 50 0) ,BL( 1 E00 ) ,B8 ( 1500 > , BR ( 1500) , CT( 1500 ) , CLt 1500 ) , CB( 1500) , 

1 CR( 1500 ) , DXT ( 1 500 ) , DXL! 1500) , DXBC 1500 ) , DXR (1500) 

NEL=N «■ I 1 

CALL CONNECT IN. C> 

DO) 11 = 1 ,114-1 

X ( I > = '/ F 1 U 4 ( I 1 ) *■ t XF 1 D - X F 1 U ) / N 
D01 1J*1,11 

11 X( J4INM 1 + I )"X!1) 

D0 1 2 1 = 1 ,n + 1 
D0 1 £1 = 1 . s 

YU=YF 1U + ( YR 1 U- YF 1 U ) •* ( J~ 1 ) X5 . 
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YD**YF 1D+(YR1B~YF1D)*( J - 1 )/S. 

Y((J-n + (N4l)+I) = YU-( 1-1 )* ( YU -YD )/N 
YU-YR2U+ l YF2U-YR2U ) * ( J- 1 )/5. 

Y D s Y R 2D *M YF2D-YR2D) J M J-1 )/S 

12 Y( ( J + 5)*MN4n+I )=YU'( 1-1 >*IYU-YD)/N 

D020 1 E= 1 ,14EL 

DXT ( IE) =1 . 0 V 

DXL { I E > *1 .0 
DXB ( l E ) - 1 .0 
DXR ( IE) =1 .0 
CT(IE)«1 0 
CHIE)«1 0 
CB( IE )-*1 .0 
20 CR ( I E ) =* 1 0 

DOE 1 IE-1 , NEL 
D022I * 1 ,4 
J**C ( I E > I ) 

E X ( I ) - X (.1 ) 

22 E Y ( I ) ~ Y ( J ) 

BTIIE ) -EY ( 4 ) -E Y 1 1 ) 

BB( IE )*EY<3)-EY(2) 

B3-EX 1 3 ) -EX I 4 ) 

B4=EYt 1 J-EY<c ) 

B5»EY ( 3 ) -EY ( 4 ) 

B6*B4+B5 
AC 1 ^BT ( IF > +B3 
AC2-1 . 0/3 . 0*tB3*B6 
AC ( IE J s’ A C* 1 + A C 2 

BL < IE )- SORT t B 3 * + 2 . 0 * B 4 * * 2 . 0) 

BR ( IE I “ BQR T ( P3 * + 2 . 0 + B5 * + 2 . 0 ) 

C IFIIH EQ 1 OR. IE EQ.N+1.0R.IE EQ . 2*N+ 1 . OR . I E . EQ .‘3*N + 1 . OR . 

C 1 IE.EQ.4+N+1 OR 1E.FQ.5 + NM OR . I E . EQ . &*N + 1 . OR . I E . EQ . N . OR . I E . EQ . 

C 1 2 * N OR IE E 9 3 . OR IE EC? 4*N OR . I E . EQ . 5* N . OR . I E . EQ . 6*N . OR . 

C 1 1E.EO 7H4JC07021 

Y 7 *= C E Y t *1 )*EY( 1 ) ) /2 0 
YB=( E Y( 3 ) + EY< 2 ) ) /g 0 
YC=<Y7+TBJ/3 0 
XC«(EX( 2) +EX< 1 7 ) SZ 0 



I F 1 IE .EC? 1 

OR 

I E 

EQ . 

N + t .OR 

IE . EQ 

2+N4 1 OR . IE . EQ . 3*N+1 . OR . 

1 

IE . EQ . 4 K4-f 1 

OR 

I E 

. EQ 

5 ♦ N 4 1 

OR 

JE 

EQ . G*N* 1 . OR . IE . EQ. 1 0*N+1 . OR 

t 

1 1 . EQ r 4 ro 1 

OR 

. ir 

. r o 

8 ♦ N ♦ 1 

cm , 

. ir 

. TO . 9 1 7 C OT 05 


J*=C< I C- 1 , 1 J 
EX1*X(J) 

E Y 1 «Y t J ) 

J = C( IE-1 , 4 ) 

EX8=X ( J 7 
EYB«Y ( J 7 

TTT»( EY8+EY J J /Z . 0 
YC1=( Y7 1+T T J/2. 0 
X C I - f E X ( I J4EXI 7/E. 0 

DH=*SQRT <<YCl-YCn*2.0*(XC1-XC)*42.0) 

CT ( IE )«(XC-XC1 T/DH 
DXT ( IE) -BH 

5 IF( IE , EG. N .OR . IE . EQ . 2*N. OR . IE E0.3*N.0R. 1E.EQ.4*N.0R. IE. EQ . 5*N 

1 . OR . IE. EQ . G*N .OR . IE . EQ. 7*N.OR . IE ,EQ. 8*N. OR. IE . EQ . 9*N. OR. 

1 IE.EQ . 1 O^N.OR . IE EQ . 1 1 *♦ N J GOT 06 
J = C( I E f 1 / 2 ) 

EX4«X( J ) 

ET4^Y ( J ) 

J-CI IEM . 3 ) 

EXS-X(J) 

EY5*T( J ) 

YBl«tEY44ET5) /2 . 0 
TCI YB+YB1 7 /Z 0 
xcisiBXtj)4cx^)/e.o 

DH- SORT U XC1 -XC ) *♦ +2 . 0 + « TCI -YC )♦=♦£ . 0 I 



ana 


CBC IE)»(XC1 -XO/DH 
DXEU IE J -DH 

6 IF? X E . GT . N ) THEN 

J = C( IE-N, 1 } 

EXH*X( J > 

EY2«YIJ) 

J-C ( IE-N, 2) 

EX3«X ( J ) 

EY3*»Y C J 1 

TT1 *t PTC + CY M ) ) y 2T . 0 
TCI - ( CY.UET < 2T y J y? . 0 
TCI »< Y7 1 + YB1 1/^.0 
DXLI XE >sr< YC-YC1 J 
CL( IE )*P3yBL( 1EJ 
END 1 F 

irt IE .LT 1 0*N)7HCN 
J«C< IC+N, 3) 

EX$*=X ( J ) 

EY6«Y ( J ) 

J*=CC IE+N, 4 ) 

EX7-X i J ) 

EY7»Y < J l 

YT 1 »( ETT + E Yt 4 ) 1/E. 0 
YB1»<EY0 + EY<3 J )/2. 0 
YC1*lY71*YB1>/£.0 
dxr ( iei-tci-yc 

CR C IE)" E3/BR ( IE ) 

END I F 

El CONTINUE 

RETURN 
END 

SUBROUTINE TO CALCULATE r HE DISCHARGE BETWEEN* THE 
NEXT TO RIVER CEILS AND THE RIVER CELLS BY USING 
THE WEIR FORMULA 

SUBROUT INF WMRUNKI 1 , K , J , I E , D X , CO , H 1 , HE , HB , Q > 

I MF LICIT DOUBLE FREC 1 S I ON ( A-H , O-Z ) 

DIMENSION DX( 1 500 ) , COC 1500) ,H!( 1 00, 1 1 ) , HH U 00, 1 ! ) 
C1®£ . / 3 . 1-SQR7 (2 0*9.81 MDX(IE) 

C2«C1 *1 . S 

IFCHSU.K) LT. HE) THEN 
U«0 . 0 
H sH8C I * J 1 
CDO = 0 . 0 
CDU=1 . 0G 
D = 1 .0 
ELSE 

IF ( ( H2 1 I » K ) -HB ) . LE . He ( I , J ) ) THEN 
W=HS< I , K 1 -*HB 
H«HH( l , J >«U 
D*= I . 0 
ELSE 

W*HSU,J) 

H*H2( I , KJ-HB-W 
D= ~ 1 . 0 
END I F 

I F l W L T < 0 867* t W*H> ) ) T HEN 
CPQ- 0 0 

CDW-0 .MHO, 0 TS'* ( H/U ) 

ELSE 

CDO" 0 6 I I -0 . 1 75 U U/( H + W ) ) 

CDW-0 . 0 
END I F 
ENDIF 

IFtH.LE . 0 . 0 ) THEN 
Q*0 . 0 
ELSE 
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Q*D*<C1 *CDU*H+t 1 .5 + C?+CDO'*U>»H++0.5) + CO( IE) 

Q1*C1 +CPU4H++ 1 $ 

Q2“C2 - »CDO + U + H'* + 0 . 5 
Q“D+ ( Q 1 + 0 2 ) + C 0 ( IE ) 

E ND I F 

RETURN 

END 

SUBROUTINE CALCULATES THE DISCHARGE COMING 
IN OR GOING OUT FROM/TO TOP AND BOTTOM CELLS 
SUBROUTINE TOPBOTIK I , J, IE,FN,RH,Z1 , Z2 , H2 , HB , BT , BB , 

1 CT,CB,DXT,DXB,GS,QI, ON , OM , QT , QB ) 

IMPLICIT DOUBLE PRE C I S I ON ( A -H , O-Z ) 

DIMENSION Z1 (100,1 1 ) ,13(100, 11 ) ,H3(100, 11 ),BT(1500),BB(1500) 
1 , CT ( 1500) ,CBt 150 0) ,DXT ( 150 0) ,DXB( 1500) 

IF( I .EG . 1 ) GOT 0 1 0 

HT*(HE( I . J )+M£( 1-1 . J ) )/£ 0 

I F< Z2( I , J ) . L T . c 2( I -1 , J 1)1 HEN 

$T»SQRI ( f Z2( I- 1 , J )-Z£( I , J ) J/DXT ( IE ) ) 

ELSE 

ST = -SQRT( (Z2( I , J1-Z2U-1 , J ) )/DXT( IE) ) 

END IF 

0 HBT-IH2U , J ) tH3( I + 1 , I ) )/2 . 0 

IF(Z2(1,J).LT .22(1+1 , J ) ) T HEN 
SB*SGRT t ( Z2( 1+1 , J ) - Z 2 ( I ,J) )/DXB ( IE ) ) 

ELSE 

SB=-5QRT( (Z2U,J)-ZeU+l,J) )/DXB( IEJ ) 

ENDIF 

I F ( J . EG . 1 . OR J EG ,11) THEN 
IF( I . EG . 1 ) THEN 
QT = Q5 
EL5E 

GT=BT< IEM*(S / 3 ) +C T ( IE M'HT * ► ( 5 . /3 . ) *ST 
1 /CFN*CBT( IEI+HD '*12 ,/3. ) ) 

ENDIF 

QB*BB ( lEM'MS, / 3 . ) * C E ( tE)*HBHM5./3. ) 4 SB 
1 /(FN'MBPl IE > *m?1 M 4 (2 .'3 >1 
ENDIF 

IF( J . GE . £ . AND J.1.E E.OF.J.GE T . AH D . J . LE . 1 0 ) T HEN 
I F ( I .EG. 1 ) THEN 

IF( J . GB £ . AND . J . L E . 4 . OR . J . GE . S . AMD . J . LE . 1 0 >GT*ai 

IF( J . EG .S . OR . J . EG . 7 ) GT»GN 

ELSE 

GT*BTC lEM'HTM (5. /3. M 5TI CTt IEJ/FN 
ENDIF 

GB = BB( IE) *HB T t e ( S . /3 . ) * SB*CB ( I E ) /FN 
ENDIF 

IF( J.EG .0 ) T HEN 
I F ( I . EG . 1 ) THEN 
at = GT1 
ELSE 

l F t HT . LT . HB ) T HEN 
P T = B T ( IF ) 42. 04HT 
ELSE 

PT*BT ( IE) +2 . OTHB 
ENDIF 

QT* ( B T ( tEM.HI)H (S./3. )*ST+CT( IE)/CRN*PT**(2./3. >) 

ENDIF 

IF( HBT .LT .HB JTHEN 
PB*BB( IF ) »c . 0-t HBT 
ELSE 

PB*BB( IE) +E . 0 1 HB 
ENDIF 

QB* ( BB t IE) < HE T ) * 4 (B . /3 . MSSiCBt IE ) /" l RN * : PB* 4= ( E . /3 . ) ) 

END IF 
RET URN 
END 
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•SUBROUTINE TO CALCULATE THE DISCHARGE COMING IN 
OR GOING OUT FROM/TO LEFT AND RIGHT CELLS 
SUBROUT INE LTRT ( I . J . K , IE . FN, Z1 . ZE, HE, BO, CO, DX,Q ) 

IMPLICIT DOUBLE PREC 1 S I ON < fl-H , O-Z ) 

DIMENSION Z1 (100,11 ),Z£ (100, 11), HE (100, 11), BO (1500) 

* , CO( 1 50 0 > , PX ( 1 500 ) 

H«(Hc( I , J)+H£( I , K> )/E .0 l 

IF t Zc( I , J ) . LT . 2E< I ,K ) ) THEN 
SL=SQRT ( 1 ZS l I , K ) - Z 2 ( I ,J) ) / D X ( IE) ) 

ELSE 

SL*-5CmTt(Z2(l,J)-Z5CI,K>)/PXME)) 

ENDir 

Q»B0tIF)*H«M5 /3. 1+5L + C0MEJ/FN 

RETURN 

END 

SUBROUTINE TO CALCULATE THE DIFFERENTIAL TERMS COTRI BUT I NG 
TO MAIN CELL FROM LEFT/RIGHT THROUGH WEIR LINK . 

SUBROUTINE D I F FUE I RL I NK ( I , K , J , I E , D X , CO , H 1 , HE , HB , DOR , DQF ) 

IMPLICIT DOUPLF PRECI SI 0N( A-H , O-Z ) 

DIMENSION DXM 509 ) .HI (100, 11 ) , H2( 1 00, 11 ) , CO( 1 500 ) 

C1*E./3.*S0RT(E. 0*9. SI > *DX ( IE ) 

CE=C1*1 5 

IF ( HE( I , K ) LT . HP IT HEN 
H e HE ( I , J) 

CDUM 0 6 
0CJR-O . 0 

DGF = C 1 * COW* I 5tHH(i,5iC0UE) 

ELSE 

I F( ( HEt J , K ) -MB ) . LE HE ( 1 , J ) ) THEN ' 

U*HE ( I , K) -HB 

H=HE( I , J ) -HEt I , K 1 H1B 

I F ( H . LE .0.0) THEN 

DGR = 0 . 0 

DQIJ” 0 . 0 

ELSE 

A1--CH0 611*1 5HHH0.5I 

A E * - C 1 *0 075+ ( l)*E . S+CM* tl 5) + ( H + + E .5) ) / ( LM< + 2 . 0) 

A3-C5+0 . 5 1 1 + (H» * 0 . S-U/(E . O+SORT (H) ) ) 

A4=-CE< 0 . 1 75UJ< ( (H+*0 ,5)+E .O-U+O . 5/SQRT(H) )/(U+H) 

A5=C1 *0 . Q75*£. 5+H++1 5/U 
A 6 =» C c + 0 6 1 1 *U* 0 . 5/ SORT ( H ) 

A7=-C2+ 0 . 17Sm'H2.0M( ( U + H ) * 0 . 5/SGRT ( H ) -SORT ( H ) )/( (U+H)**E. 0 ) 
IF(U.Lr.(0.6S667*(U+H)> ) THEN 
DQR=(A1+A5)+CO( IE) 

DQF ®(-A1+A5)+C0( IE) 

ELSE 

D0R“( A3 + A4 ) +CO< IE ) 

DOF=(A6+A7)*CO( IE) 

END IF 
END IF 
ELSE 

U*H5( I , J ) 

H=H2 ( I , K ) -HB-H? t I , J ) 

I F ( H . LE . 0 . 0 ) THEN 
DQR- 0 . 0 
DQF =0 . 0 
ELSE 

A1=C1 *0 . 6 1 1 +1 . 5+H* <0.5 
A2 = C1 *0 .0 75 + 5 .5* (H++ 1 . 5)/U 
A3=C5< 0 .61 1 +U + 0 . 5/ SORT (HI 

A4=-CE'0 . 17S+(U+*2 . 0 )*( (U+H) *0 . 5/SQRT ( H ) -SORT ( H ) ) / ( ( U+H ) **£ . 0 ) 
A5=-C1 *0 . 075+ (U+E .5+ CH.+ *1 . 5 ) + ( H+ +£ . 5 ) ) /( U**E . 0 ) 

A6=C£+0 .611 + ( H * * 0 .5-W+0 .5/SQRTCH) ) 

A7-»-C£ < 0 1 75* U + (£ . Q+SORTC H )-U+0 . 5/SQRTCH) )/(U+H) 

I F ( U LI . ( 0 . 66667+ (U+H) )) THEN 
DOR * < -A 1 -A£MCOC IF. > 


ne* 



DQF*( At -AS ) *C0( IE ) 

ELSE 

DQR = (-A3-A4> * CO C IE) 

DQF = (-A6-A7HCO( IE) 

ENDIF 
END IF 
ENDIF 

END IF » 

RETURN 

END 

SUBROUTINE TO CALCULATE THE DIFFERENTIAL TERMS OF DISCHARGE 
CONTRIBUTING TO THE MAIN CELLS FROM LEFT AND RIGHT CELLS . 

SUBROUT INF D I FFL T R T C I , J , K , IE.FN.Z1 , Z2 , HE , BO , CO , DX , DQ 1 , DQE ) 

IMPLICIT DOUBLE PRE C I SI ON t A -H , O-Z ) 

DIMENSION Zt ( t 00 . t 1 ) , Z2( 1 00, 1 1 ) , H2 ( 1 00, 1 1 ) , BO( 1 500 ) 

* , COU500) ,DX( tSOO) 

H*H2( I . J )+H2( I ,K > 

IF(Z2U,J).LT 72(1 , K ) IT HEN 
SL-SORT (?2< I , K J-7£('I . J ) ) 

DSL 1*--t .0/(2 0 *-SL ) 

D5L2® I . 0/(2. 0 *SL ) 

ELSE 

SL--SORT ( 72 ( I , J > Z2< I ,K) 1 
DSL 1 - 1 . 0/T2.0»SL ! 

DSLS' 5 - 1 .0/(2. 0«SL) 

ENDIF 

A1 =BD( IE )/(2. 0 M- (5. /3 . ) ♦•FN t'SQR T ( DX ( IE) ) ) 

I F ( AB5 ( 5L ) . LE .0.0 ) THEN 
D01 =0 . 0 , 

DQc = 0 . 0 
ELSE 

D0 1 * A 1 * ( SL M 5 / 3 . M H » ♦' ( 2 . / 3 . )tHM(S,/3. M'DSLI ) 

DO£»A1 «' ( SL *• (S . / 3 )*HN’(E./'3.)*tUt(5./3. ) *DSL2 ) 

ENDIF 

RETURN 

END 

SUBROUTINE CALCULATES THE DIFFERENTIAL ELEMENT5 OF 
TOP AND BOTTOM CELLS CONTRIBUTING THE MAIN CELL . 

SUBROUt INF. OIFf'TOPBO TM( l , J , IE, FN,RH , Z1 , ZZ , HE, HB , BT , BB , CT , CB , 

1 DXT,DXB,PQT1<DQT2,DQB1 , DGB2 ) 

IMPLICIT DOUBLE PREC I S I ON ( A-H, O-Z ) 

DIMENSION Z1 ( 1 00 , 11), Z2 (100, I 1 J , H2 (100, 11) 

* , BT (1 500) . CTU500 ), BB( 1500 ) , 

t CBN 1500) ,DXT( 1500) ,DXB< 1500 ) 

IF( I . EQ . 1 )G0T05 
Hf»))2( I . .1 ) Mir ( 1-1 , I) 

1 F ( Z? ( I . J ) L I . /2 ( 1 • l < J • ) 1 1lf N 
ST*SQR T ( 721 1-1 , J )-Z2 ( I , J ) > 

DST 1 * - 1 0/(2. 0 +ST ) 

DSTH-l 0/ ( 2 . 0 <• ST ) 

ELSE 

ST s =-SQRT(Z2( I , J)-Z2( I - 1 d) ) 

DST1*1 .0/(2 . 0< ST ) 

DST ?=- 1 0/(2. 0+ST) 

ENDIF 

HBT=H2( I . J)+H2(I+ 1 . J > 

I F( Z2( 1 , J ) • L T • 1 + 1 . J > ) THEN 

SB»SORT( ?g( 1 + 1 ,J)-i£U »JI ) 

DSB 1 *- 1 . 0/(2 0*5P) 

DSB2— 1 • 0/(2. 0+ SB ) 

ELSE 

5B--50R T ( Z2 ( 1 . J ) - Z2 ( I + I , J J ) 

DSB1 =1 . 0/ < 2 0 ) 

03BS— - T 0/(E.0'5B) 

END IF 

I f t J EO 1 . OR J .20 . T 1 ) T HEM 
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IF( I . EQ 1 ) THEN 
DQT 1*0 . 0 

dote- o . o 

ELSE 

A 1 * 8T I IE > *«' t S . / 3 M Cl t IE )/l FN*E . O+SORTl DXTl IE) ) > 

A3*2 0 + BUIE) + HT 

AS*t5TM5 / 3 MHMUS./J ) + H T * * C 5 . /3 . ) TOST 1 )/A3**t2./3. ) 


A6*HT ♦ * r ? / 3 )*ST*£.0/13 0 t A3M( 5 . /3 . ) ) 

AS* t 5 T * t S / 3 J*Hr*M£ /3 . ) + HT» ♦ l 5 . / 3 . ) *DST2 ) / A3 ** C 2 . /3 . ) 
DOT 1»A) ♦ ! A5-A9 I 


DQTt*A) * C AS-AG 1 
END IF 

A2*8B( IEMM5./3 IKBl IE ) / ! FN18 . 0*3CIRT (DXBt IE) )) 

A-<) = 2. 0 1 BB t IE1 + HBT 

A7«t5B*t5./3 )*HBTT*t2./3. >+HBTt*t5./3. ) *DSB1 ) /A4 * * l 2 . /3 . ) 

A8-HBT*4I5./3. JtSB^B 0 / I 3 . 0 * A4 %* t 5 . /3 . ) ) 

A10»ISB'*>l5./3 MHBTTTIE./3. >+HBT**l5./3. ) *DSB£ ) /A4* * t 2 . /S . ) 
DOB I *AE* I A7 - A8 ) 

DQBS*AE*t AI 0-AS) 

END I F 

IFJ J . GE . ? . AND . J . LE . 5 DR . J . GE . T . AND . J . LE . 1 0 ) THEN 
A1*BT.t IE) 4-CT t IE ) / ( FN *2 0* * IS ./3. JTSGRTtDXTl IE) ) ) 

AE = BB t IEMCSf I E ) / t FN«E 0 * * t 5 . S3 . ) *50RTC DXBt IE) ) ) 

IF! I .EQ 1 ) THEN 
DQT 1*0 0 
DQTE*0 . 0 
ELSE 

DQT 1 *A1 *’ I ST* ( 5 . /3 . ) *HT ♦ * I 5 . /3 . ) 4 H T* * t 5 . /3 . ) *DSTI ) 

DQT2=A1 * t STf 15 . / 3 . ) *HT**l E . / 3 . ) + H T* * t S . /"3 . ) *DST2 ) 

ENOIF 

DQBT a A2'MS0M5 /3. MHBT* Me. / I . ) + HB T ** I 5 . /3 . ) *DS8 1 ) 
DaB£*A£*-t5B l M5./3. )*HBT*M2./3. ) H)BT + * (5 . /3 . ) *D5BZ ) 


ENDIF 

IFt J .EQ . $ ) THEN 
I Ft I .EQ . 1 >THEN 
DQT1 *0 . I) 

DQT2*0 . 0 
ELSE 

A1-BT l !E> ♦* tS . /3 . > *CT l IE)/f RN* E • 0 M t5 ./3. ) *SQRT t DXT ( IE) ) ) 
A3*ST*MS./3. ) ♦HTt * ( 2 ./3 . )+HT+-t- (5. / 3. ) *DST1 
A<4*ST*(5./3. >*'HT< *(2 / 3 . )+MT *♦ (S./3. )*DST2 
I Ft (HT/2 0) .LT.HB) THEN 
PT*BI ( IE) +HT 

AS* ( HT *' * f 5 . / 3 )*STM2./3 ) )/PT*» ( S . /3 . ) 

A6»A3/PT*+t2./3. ) 

A7*A4/PTff(£./3. ) 

DQT1 »A1 * ( A6-AS ) 

DQT2-A1 *( AT-A5 ) 


ELSE 

PT»BT( IE>+2. OiHB 
A5°A3/P T *’ * ( £ ./3. ) 
A6*A4/PT <‘*t 2 . /3 . ) 


DQT1 =A1 *AS 
DQT2*A1 »A6 


ENDIF 

A2 = BB t IE) * MS . / 3 ) « C E t IE ) / t PN< £ 0tHS./3. XSQRTtDXBt IE) ) ) 
A8*SB* - t c . / 3 ) * HB T ♦ i ( 2 . / 3 ) *HBT » * C 5 . /3 . ) *DSBt 
AS=SB*<5 / 3 1 )*HBT*M2./3 ) -»HBT * * ( S . /3 . ) » DSB2 


I F t (HBT /<? . 0 ) LT.HB) THEN 
PB-BE( IEJ+HBT 

A5=HBT * *15 / J. )*SB) 12./3. >/PBM (5./3. ) 


A6=AB/PB) t(£./3.) 
A7=*AQ/PB» H2./5 ■ ‘ 


DOB 1 *AS + 1 A6-A5 ) 


DQB2»A2’M A .’-AS) 
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5 


ELSE 

P8=8B( IP.) +2 . O + HB 
AS-A8/PBM ( £ /3. > 
A6*A9/F B + M2 / 3 ) 
DQBt “AS ♦ AS 
0QB£-»A£+A6 


ENDIF 
END I F 
RETURN 
END 

ELEMENTS OF THE JACOBIAN MATRIX 

SUBROUT INE J ACOB I N . JN , FN , RN , Z1 , Z2 , HI , HE, DT, AC, OXT , DXL, DXB, 

1 DXR , BT , BL ,BB , BR ,CT , CL ,CB, CR , MB , OF ) 

IMPLICIT DOUBLE PREC I S I ON ( A-H , O-Z ) 

D1MET4S1 ON 21(190, 1 1) ,H1 I 1 00, ITT, HE II 00 , 1 1 ) , ACI 1 50 0 ) , DXT 1 1500), 
1 DXL ( ISO 0 ) ,DX8t 1 SOO ) , DXR < 1 SOO ) , BTI 1500) ,BLM500),BB11500), 

I CT { 1 SOP ) , CL f I SOO ) , CBI 1 SOO ) , CR ( 1 50 0 ) , DF ( JN , JN ) , 22 ( 1 00 , 11 ) 

* , BRMSQO) 

00251-1.14-1 

IE-™I 

CALL DIFF TOPCOTMl IE, T , IE. FN, RN, Z1 , ZE , HE , MB , BT , BB , CT , CB , DXT , 

1 DXB , DOT 1 , DOTS , BQB1 , DOBS ) 

CALL DIFF l TRT I IE . 1 ,E, IE ,FN , Z1 , ZE, HE, BR , CR, DXR, D01 , DQ2) 

D0E5 J = 1 , JH 

I FI J ED I OR J EG). I + N-1 .OR . ( I .LT.N-1 .AND. J.EQ. 1 + 1 ) -OR. 

1 II.GT l.AND.J ED 1-1)) THEN 

IFCJ ED I )DF I I , J ) * ACI IE ) /DT-DDB1-DQT 1 -DQ1 
IF! J . EO . I + 1 )DF< I , J )=-DQB2 
IF! J ED . 1-1 )DF( I , J )=-DOT£ 

IFIJ EO I *-N-l )I»F( I , J ) k -D0S 
ELSE 

DF I I . J) -0 0 
ENDIF 
CONT INUE 


1 


1 


JO J»E, <1 

30261 *J ♦ ( N- 1 1-N+S, JMN- 1 > 

JE-I+J-l 

CALL D1 FP TOP POT MI I 1 , J , IE , FN, RN , Z1 , ZE , HE , MB , BT , BB , CT , CB , DXT , 

DXB, DOT 1 , DOT E , DOB 1 ,DOBE) 

C^LL DlFPLIRnn,J.JM,IE,FN,Z1,Z 2 ,H 2 ,BR,CR,DXR,DaR1,DQR2 

CALL DIFF l TRT I II , J , J-T , IE, FN, Z1 ,ZE,HE, BL, CL, DXL, DQL1 , DDLS 


I FIJI .ED. I OR.J1 EO.l-W-1.OR J1.E0 I -N+ 1 . OR . I I . LT . J* I N- 1 ) . AND . 
J1 ED I + 1 J . OR . I I . G T . J * I N-1 )-N + E . AND . J1 . ED . 1-1 ) ) THEN 
IF! J1 .ED. I )DF( I, J1 )=ACI IE ) /DT-DOB1 -DDT1 -DDL1 -DQR1 


IFIJI .FO. I + 1 ) OF I I , J 1 )— D0B2 
I F I J V . ED . I- 1 > OF ( I , J 1 ) --DOTS 
IFIJI .EO. 1 ♦N-1 * DF I I ,J1 ) =-DOR£ 
IFIJI . EO I“N+ 1 )DF( I, J1 )=-DQL£ 


ELSE 

DFI I , J1 )*0 0 
ENDIF 
CONTINUE 
ENDDO 

DO£T I '*•1 *N-S, *" * N - 5 

I 


1 


1 


LL + Dirrior BO 1 M < I 1 , 5 , I E ■ F N , RN , Z T , zs , HE , HP , BT , BB I CT , CB , DXT , 
LL°DI F^l'tRT < n B S ', 4? IE ! FM . Z 1 , ZE , HE , BL , CL , DXL , DGL 1 , DOLE ) 

LL DIFFWEIRLINM(I1,6.S.IE.BR,CR,H1,HE,HB,DQ1,DQE) 


DOST J®! ,JN 

IF I J ED I OR J ED. I - N ♦ T • 
j c<? 1-1 ) .OR II LT s,n "' s 
IKl J ED t ) OF I I • J)-ACt IE) 


OR J ED . I + N-1 . OR . I I • GT. 4*N-3 . AND . 

AND J . SO. 1+ 1 ) ) THEN 
.-UT-ODBT -DDTT-0DL1 +DDE 
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IFtJ EG I + 1 ) OF ( 1 , J 1 *> -DOB? 
IFtJ EG 1-1 M>F< I , J )*-OQT£ 

1 F( J . EG . UN- I )OFU , J )=DOI 
IFtJ. EG I-N+1 >OF l I , J )»-DQl2 
ELSE 

DFU , J 1 -0 0 
END I F 


1 


1 


1 


1 


C ON T INL'E 


D02SI>-5»N~4 f 

1 1*1 -5*H*5 
1 E* I +5 


CALL PlPPTOPBdni II ,0, IB , FN , RN , Z1 , ZZ , HZ , HB , BT , BB , CT , CB , DXT , 
DXB , DOT! . DQTS , DOB 1 , DOBS ) 

CALL DIFFUEIRt. INK Ml , 6 , 5 , I E , EL , CL , H 1 , HE , HB , DGL 1 , DOLE J 
CALL DIFFUEIRLINKl II ,6,7, IE,BR,CR,H1 , H2,HB,0QR1 , DQRE ) 

DOE8J* J , JN 

IFtJ. EG. I OR J EG. I-N+1. OR J .EG. I+N-1 ,0R. I I . GT . 5+N-4 . AND . 

J.EIJ I-1> OR U L T . 6 +N- 6 AMD J . EQ . I + 1 ) ) THEN 

IF t J EIJ.DDFU , JHACUEJ/DT-OQB1 -DQT1 -DGL1 -DQR1 

1FI J .EG . J +1 IDFl I , J)»-DGB2 

IFtJ EG 1-11DFU f J)*-DGT2 

1 F ( J EG I 01-1 >OFU ,J )*-DQR2 

IFtJ EG I -M+1 >DFU , J )*-DQL2 

ELSE 

DFt I , J ) 0 

ENDIF 

CONTINUE 

Doeei*s+M-s, r + N- r 
! I»!-6MR{ 


ir.*i +t? 

CAL L 01 FT TOPpOtrU I 1 , 7, IE , FN ,RN ,Z1 , Z2 , HS, HB, BT , BB , CT , CB , DXT , 
DXB , DOT 1 . OGT 2 , DOB 1 , DOBS) 

CALL OIFFLTR t < H . T , 8 , 1 E , FN , 7 1 , 22 , H2 , BR , CR , DXR , DQR1 , DQR2) 
CALL DIFFUEIRLINIt t 1 1 , 6 . 7 , IF. , BL , CL , H 1 , HB , HB , DQ1 , DUE ) 

D02?J*1 , JN 

IFtJ EG l OR J EG. I-N+1. OR J F.O I + N-1 . OR . ( I .GT . 6*N-5 . AND . 

J EG I-ll.OR (I.LT 7 *N-7 . AND . J .EG . 1 + 1 ) > THEN 
IFtJ EG I >OFt I . J1*ACUE)/0T-0GB1-DQT1-0QR1+DQ£ 

IFl J FO, 1+1 1 OF ( 1 , J)*-DGB2 
IFtJ. EG 1-1 HJFU . J)=-DQTE 
IFtJ.FG. ltll'l J DF < I,J)*-D9R? 

I F ( J . EG . I-N+1 >DF< I , J )=DQ1 


ELSE 

DFt I , J 1*0 0 
ENDIF 
CONTINUE 
DO J*8. 1 0 

002 3 1 - J * t N- 1 ) - N + 2 , J * ( N- 1 ) 


IE* 1+ J-1 

, „„ 

CALL D IF F TOP BO T 11 ( 1 1 , J, IE .FN , RN, 21 , Z2, H2, HB, BT, BB, CT, CB , DXT 

?^ D S!;n:?Ruri:^M! E . f N,Z-, 2 ,,M S ,BR,CR,DXR,D 0B ,,D O « S . 
CALL DIFFLTRTtll . J , J — T < IE , FN, Z1 , Z2 , H2 , BL , CL ,DXL/ DQL1 , DQL2) 

IFtjrto’rOR J1 60 I4N-1.0R.J1 EG. I-N+1 . OR . ( I . LT . J* ( N- 1) . 

1 !,! OR U 61.JKM-1WH2 AND.J1.E0.I~1I) THEN 

IFtJ! ! F 0 I»om, J1I-ACUE)/DT-D0BI-D0T1-DQL1-DQR1 
IFtJI EG I + 1)DF(I,J1)*-DGB£: 

IFtJI FG I - 1 IDF ( I , JM )*-D0T2 
IFtJ! F.O I+N- IlDFt I, J1 >*-D9R8 
IFtJI. F 0 I * N * I >Br- 1 1 , J I ) a "DGl 2 


ELSE 

om.jD’O.o 

ENDIF 

CONTINUE 


AND. 
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24 


ENOOO 

D024I - 1 0 *N-3 ,11 + <N-1 ) 

II * I - 1 O’* N+ 1 0 
IE* 1+ 1 0 

CALL D1FF TOPBOTtlt 1 1 , 1 1 . I E , F N , RN , Z t , ZS , HS , HB , BT , BB , CT, CB, DXT, 
1 DXB , DOT 1 , DOTS , DGB1 ,DQB2) 

CALL DIFFLTRT 1 II ,11 , 1 0 , I E , FN . Z 1 , ZE , H£ , BL , CL , DXL , DQ 1 , DQ2 ) 
D024J»! , JN 

1F1J EQ I OR J EO.I-N+1 OR (I LT . 11 *N-1 1 . AND . J . EQ . 1 + 1 ) . OR . 

1 tl GT.10«N-9 AND J EQ. 1-1)) THEN 

IF 1 J EO . DDF ( l „ J >=AC< IE) /D1 -DQB1 -DQT1 -DQ1 

IF(J FQ I +1 )l»F < I , J )*-DQB2 

IF1J.F.Q, 1-1 )OFt I ,J)*-DQTE 

1 F ( J EQ . HIM >DF 1 I , J )*-DQ2 

ELSE 

DF ( I . J ) ~0 0 
END IF 
COM Tl HUE 
RETURN 
END 



